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1.  Overview 


Equilibration  of  entangled  polymers  is  nontrivial  even  for  a  case  of  coarse  grained  chains 
because  of  slow  reptation  dynamics  exhibited  by  high  molecular  weight  chains,  where  center-of- 
mass  diffusivity,  D  scaled  with  polymer  length,  N  as  D  oc  AT2'3 .  Even  with  modem  day  parallel 
computers,  brute-force  equilibration  of  well-entangled  branched  polymers  is  prohibitively 
expensive.  Standard  method  in  which  one  starts  with  an  ensemble  of  chains  with  the  correct  end- 
to-end  distance  arranged  randomly  in  the  simulation  cell  and  introduces  the  excluded  volume 
rapidly,  leads  to  deformation  on  the  intermediate  length  scales.  Another  set  of  long-standing 
methods  includes  the  chain  connectivity  altering  algorithms.  These  bridging  methods  employ  a 
set  of  complex  Monte  Carlo  (MC)  moves  that  can  be  difficult  to  implement  for  branched 
polymer  architectures  and  these  MC  moves  have  very  low  acceptance  ratio  (i,  2). 

New  preparation  protocol  has  been  established  to  overcome  the  described  difficulties.  This 
protocol  is  general  to  any  polymer  architecture  and  does  not  bring  any  deformation  of  polymer 
chains.  Short  simulation  runs  with  soft  potential  used  in  Dissipative  Particle  Dynamics  (DPD)  (3) 
is  performed  for  an  ensemble  of  polymer  chains  with  the  correct  end-to-end  distance.  After  a 
gradual  increase  of  the  strength  of  DPD  potential,  short  simulation  with  target  coarse-grained 
potential  is  performed. 

In  this  report,  we  adopt  our  method  for  two  coarse-grained  models  1)  Kremer-Grest  model  ( 4 ) 
and  2)  DPD  model  with  Segmental  Repulsive  Potential  (mSRP)  (5).  However,  the  described 
protocol  is  general  and  applicable  for  any  coarse-grained  model. 

The  present  preparation  method  exploits  two  programs:  gel  and  Large-scale  Atomic/Molecular 
Massively  Parallel  Simulator  (LAMMPS).  Gel  is  an  in  house  written  C  program  of  coarse¬ 
grained  polymer  builder,  and  LAMMPS  is  a  molecular  dynamics  program  from  Sandia  National 
Laboratories  (6). 


2.  Models 


2.1  Kremer-Grest  Model 


Kremer-Grest  model:  The  pair  interaction  between  topologically  nonconnected  particles  is 
described  by  the  standard  truncated  Lennard-Jones  pair  potential: 


Uu(r)  =  Ae 
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,  where  s  is  the  depth  of  the  potential  well  and  a  is 


the  distance  where  interparticle  force  is  zero.  rc  represents  the  cutoff  distance.  rc  =  21/6<x  is 
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chosen,  yielding  so-called  Weeks-Chandler-Andersen  excluded  volume  potential,  Uwca- 
Topologically  bound  monomers  interact  according  to  the  standard  FENE/Lennard- Jones  bonded 
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,  where 


standard  parameter  values  for  of  Ro  =  1.5  a  and  a  =  30  s/a2. 


To  model  different  flexibility  of  polymer  chains,  a  bending  potential  of 

form  Ubend(r)  =  kd  (l  -cos  6)  is  employed,  where  6  is  the  angle  between  two  neighboring  bonds 

and  ke  is  the  bending  constant.  The  detailed  model  description  can  be  found  in  the  original  paper 
{4). 

DPD  model  used  as  the  first  step.  The  pair  interaction  between  topologically  nonconnected 
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particles  is  described  by  soft  potential  UDPD (r)  =  —  1 - ,  where  aBpd  =  25  kBT and  rc=  1  in 

2  l  rc) 

DPD  units  of  length  are  standard  DPD  parameters  (5).  Topologically  bound  monomers  interact 
according  to  the  standard  FENE/Lennard-Jones  bonded  potential,  Ufene/lj- 
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,  where  standard  parameter 


values  for  of  Ro  =  1.5  a  and  a  =  30  s/a2. 

2.2  DPD  Polymer  Model  with  mSRP 

The  pair  interaction  between  topologically  nonconnected  particles  is  described  by  soft  potential 


U  DPD  (r)  — 
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where  adpd  is  the  maximum  repulsion  in  kpT  units  and  rc=  1  in  DPD 
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units  of  length.  Topologically  bound  monomers  interact  according  to  the  harmonic  potential, 
Uharm-  U harm{r)  =  ~^{r - r(}  f ,  where  k is  the  harmonic  constant  and  ro  is  the  equilibrium  bond 

distance.  In  order  to  avoid  bond  shrinkage  and  make  a  distance  between  bonded  particles  close  to 
the  distance  between  nonbonded  particles,  the  values  of  k  =  225  kBT/p  and  ro  =  0.85  rc  are  used. 
Chains  crossing  as  a  result  of  soft  DPD  potential  is  greatly  reduced  by  introducing  the  Segmental 
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Repulsive  Potential  (mSRP),  USRP(r)  =  — ^  1 - .  The  potential  is  effective  for  any  pair  of 

2  l  dc) 

bonds  of  the  distance  between  midpoints  of  interacting  bonds,  d  within  cutoff,  dc  with  exception 
of  neighboring  bonds  (no  repulsion).  The  repulsion  is  equally  divided  between  beads  composing 
the  bond.  The  values  of  asRp  =100  kBT and  dc=  0.8  rc  significantly  decrease  the  number  of 
topological  violations.  As  a  result  of  absence  of  bond  repulsion  between  two  neighboring  bonds 
and  soft  nonbonded  potential,  the  value  of  average  bond  angle  9  becomes  greater  than  n/2. 
Addition  of  bending  potential  Ubend{r)  =  ke{\  -  cos  6)  resolves  this  problem.  The  value  k0=  2 
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kBT  is  enough  to  compensate  the  action  of  mSRP  potential  and  polymer  behaves  as  a  flexible 
chain.  In  depth  explanation  can  found  is  in  the  original  paper  (5) 

During  all  steps  of  our  melt  preparation  we  have  used  DPD  thermostat  on  pairwise  interactions. 
In  a  DPD  simulation  of  polymers,  particles  interact  with  each  other  via  a  pairwise,  two-body, 
short-ranged  force,  F,  that  is  written  as  the  sum  of  a  conservative  force,  Fc,  dissipative  force,  F°, 
and  random  force,  F*,  as  follows:  F;  =  2X  +  +  2X-  For  a  more  detailed  description  of 

j*i  j*i  j*i 

the  DPD  thermostat,  see  the  original  papers  by  Hoogerbrugge  and  Koelman  (7)  and  Espanol  and 
Warren  (8). 


3.  Generation  of  Initial  Structure  Using  gel 


3.1  Theoretical 

To  equilibrate  a  polymer  structure,  it  is  necessary  to  generate  an  ensemble  of  polymer  chains 
with  the  correct  end-to-end  distance  (7). 

Polymer  chains  are  generated  inside  a  periodic  cubical  box  with  the  correct  end-to-end  distance 
R2(N )  =  12Cn(N—1)  on  large  length  scales  (N  is  large),  where  TV  is  a  number  of  beads,  /  is  a 
bond  distance  and  CN  is  Flory  characteristic  ratio  and  CN  =  Cx  at  TV  — >  co .  For  a  Gaussian  chain 

C  -1 

— — .  For  example,  fully  flexible 

C.+l 

Lennard-Jones  chains  with  rc  =  21/6cr  have  Cx  =1.88  and  DPD/SRP  polymers  have  Cr  =  1.32 . 

Chains  are  built  via  a  random  walk  in  3D  space  with  either  a  desired  angle  or  a  random 
angle  a  =  n  -  0 ,  that  is  sampled  around  (cos  a) .  Bead  number  density  of  Kremer-Grest  system, 

_ O  _ 'J 

p  =  0.85a  and  of  DPD/SRP  system  p  =  3rc  ,  which  is  a  standard  density  for  coarse-grained 
polymeric  melts  for  these  models,  respectively. 

3.2  Software  Description 

Program  gel  generates  melts,  blends  and  solutions.  Polymer  chains  may  include  several  blocks  of 
different  types  of  beads  and/or  angles.  Different  types  of  beads  are  used  to  drive  a  phase 
separation.  Possible  architectures  include  unconnected  beads,  linear  chains,  star-branched 
polymers  with  equal  and  unequal  arm  lengths,  combs,  dendrimers  and  hybrid  star/comb/denrimer 
structures.  Chains  are  modeled  as  flexible,  semi-flexible  or  rigid.  Figure  1  shows  the  schematic 
picture  of  polymer  architectures. 


with  a  small  bond  angle  0  the  average  cosine  is  (cos  0) 
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Figure  1.  Schematic  representation  of  bead-spring  branched 
block  copolymer.  Different  colors  correspond  to  the 
different  types  of  beads.  The  picture  corresponds  to 
example  1. 


Source  code:  gel.c. 

Arguments:  Lx,  Ly,  Lz,  R,  Nrep.  Where  Lx,  Ly,  Lz  are  the  dimensions  of  a  periodic  simulation  box, 
R  is  equilibrium  length  of  the  bond  and  is  Nrep  a  number  of  replicas.  Value  of  R  is  taken  with 
respect  to  the  bond  potential,  for  example  FENE  bond  has  R  ~  0.9655er. 
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Input:  input  is  a  topology  file  (written  by  a  user) 


Input  file 

Commentaries 

#molecules  □ 

Number  of  molecular  species 

#  atomtypes  Qt 

Number  of  bead  types 

#angletypes  □ 

Number  of  angle  types 

#polymerl  □ 

Number  of  molecules  of  polymer  1 

polymerl  □  □ 

Number  of  arms  (sidechains)  in  a  polymeric  units  at  the  same  branch  point  and  number  of 
repeated  units 

chain  1  ftp, ft';; 

Number  of  blocks,  repeated  units  and  an  indicator  of  subbranches  in  the  first  arm  (sidechain). 
The  first  arm  always  represents  a  main  chain  (backbone).  The  indicator  equal  0  means  that 
there  are  no  subbranches.  The  indicator  =1  means  that  there  are  subbranches  at  the  end  of 

□ 

this  chain  (branch  point). 

□  fj  □  □  □ 

Type  and  number  of  beads  and  cosine,  case  and  type  of  angle,  a  in  the  first  block.  Case  =  0 
is  a  random  angle  (with  average  cosine),  case  =  1  is  a  regular  angle  and  case  =2  is  a  special 
case  for  barbed  wire  (  all  3  angles  are  120°) 

branch  1  □  [jp 

Number  of  arms  (sidechains)  in  a  sub  star  and  number  of  repeated  units 

subchain  1  □  □ 

Number  of  blocks,  repeated  units  and  an  indicator  of  subbranches  in  the  first  substar.  The 

□  □□□□□ 

first  arm  always  represents  a  main  chain.  The  indicator  is  always  0. 

Type  and  number  of  beads  and  cosine,  case  and  type  of  angles  in  the  first  block.  Case  =  0  is  a 
random  angle  (with  average  cosine),  case  =  1  is  a  regular  angle  and  case  =2  is  a  special  case 
for  barbed  wire  (  all  3  angles  are  120  °) 

Output:  dataj. chain  and  snapj.xyz  files,  where  i  =  0 ,...Nrep  -1 .  data  J. chain  is  topology  file 

for  LAMMPS  software,  snapj.xyz  is  xyz  file  for  Visual  Molecular  Dynamics  (VMD)  software 
(for  a  visual  checking).  VMD  is  a  molecular  visualization  program  for  displaying,  animating,  and 
analyzing  large  biomolecular  systems  using  3-D  graphics  and  built-in  scripting. 

To  run  the  code,  a  user  has  to  have  gel.c  and  input  files  in  a  working  directory,  gcc  compiler 
could  be  used  for  compiling,  as  gcc  gel.c -o  gel  -Im.  The  user  has  to  write  input  to  define  a 
desirable  topology. 

In  the  following  example  the  code  generates  3  replicas  of  a  blend  of  linear  flexible  chains, 
composed  of  500  chains  of  250  beads  and  250  chains  of  500  beads.  The  periodic  box  size  is 
66.5x66.5x66.5  and  bond  length  is  0.9655. 

#molecules  2 
#atomtypes  1 
#angletypes  0 
Upolymerl  500 
polymerl  1  1 
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chainl  110 
1  250  -0.303135889  0  0 
#polymer2  250 
polymer2  1  1 
chainl  110 

1  500  -0.303135889  0  0 

To  ran  the  user  has  to  type  ./gel  66.5  66.5  66.5  0.9655  3.  This  action  will  produce  the  following 
output  files:  data_0. chain,  data_l .chain,  data_2. chain  and  snap _0.xyz,  snap_l.xyz,  snap_2.xyz. 

The  following  are  examples  of  input  files  for  several  possible  branched  structures: 

Example  1.  “Barbed  wire”  polymer  (see  figure  1). 

#molecules  1 
#atomtypes  2 
#angletypes  2 
#polymerl  1 
polymer  1  2  2 
chainl  110 

2  2-0.303135889  0  0 
chain2  211 

1  1  -0.303135889  0  0 
1  1  -1.01  1 
branch  1  2  1 
subchainl  210 
1  1  -0.5  2  2 
1  1  -1.01  1 
subchain2  210 

1  1  -0.5  2  2 
1  1  -1.01  1 

Example  2.  Dendrimer.  Some  of  angles  are  set  to  180°  to  make  visualization  easier  (see  figure 

2). 

#molecules  1 
#atomtypes  5 
#angletypes  2 
#polymerl  1 
polymer  1  6  1 
chainl  111 
1  10-1.0  1  1 
chain2  111 
1  10-1.0  1  1 
chain3  111 

2  5-1.0  1  1 
chain4  111 
2  5-1.0  1  1 
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chain4  111 
4  10-1.0  1  1 
chain5  111 
4  10-1.0  1  1 
branch  1  3  1 
subchainl  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4  -0.303135889  0  0 
subchain2  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4-0.303135889  0  0 
subchain3  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4-0.303135889  0  0 
branch2  3  1 
subchainl  310 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4  -0.303135889  0  0 
subchain2  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4  -0.303135889  0  0 
subchain3  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4  -0.303135889  0  0 
branch3  5  1 
subchainl  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4  -0.303135889  0  0 
subchain2  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4-0.303135889  0  0 
subchain3  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4-0.303135889  0  0 
subchain4  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
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3  4  -0.303135889  0  0 
subchain5  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4-0.303135889  0  0 
branch4  5  1 
subchainl  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4  -0.303135889  0  0 
subchain2  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4  -0.303135889  0  0 
subchain3  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4  -0.303135889  0  0 
subchain4  310 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4  -0.303135889  0  0 
subchain5  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4  -0.303135889  0  0 
branch5  2  1 
subchainl  3  10 
5  1  -0.5  1  2 
5  1-1.0  1  1 
5  4-0.303135889  0  0 
subchain2  3  10 
5  1  -0.5  1  2 
5  1-1.0  1  1 
5  4-0.303135889  0  0 
branch6  2  1 
subchainl  310 
5  1  -0.5  1  2 
5  1-1.0  1  1 
5  4  -0.303135889  0  0 
subchain2  310 
5  1  -0.5  1  2 
5  1-1.0  1  1 
5  4  -0.303135889  0  0 
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Figure  2.  Bead-spring  dendrimer.  Different  colors 
correspond  to  the  different  types  of  beads. 

The  picture  corresponds  to  example  2. 

Example  3.  Comb/dendrimer.  Some  of  angles  are  set  to  180  and  120°  to  make  visualization 
easier  (see  figure  3). 

#molecules  1 
#atomtypes  6 
#angletypes  2 
#polymerl  1 
polymer  1  710 
chainl  110 
6  40  -0.303135889  0  0 
chain2  111 
1  10-1.0  1  1 
chain3  111 
1  10-1.0  1  1 
chain4  111 
2  5-1.0  1  1 
chain5  111 

2  5-1.0  1  1 
chain6  111 
4  10-1.0  1  1 
chain7  111 
4  10-1.0  1  1 
branch2  3  1 
subchain  1  3  10 

3  1  -0.5  1  2 
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3  1  -1.01 1 
3  4  -0.303135889  0  0 
subchain2  310 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4  -0.303135889  0  0 
subchain3  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4  -0.303135889  0  0 
branch3  3  1 
subchainl  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4  -0.303135889  0  0 
subchain2  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4-0.303135889  0  0 
subchain3  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4-0.303135889  0  0 
branch4  5  1 
subchainl  310 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4  -0.303135889  0  0 
subchain2  310 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4  -0.303135889  0  0 
subchain3  310 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4  -0.303135889  0  0 
subchain4  310 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4  -0.303135889  0  0 
subchain5  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4  -0.303135889  0  0 
branch5  5  1 
subchainl  3  10 
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3  1  -0.5  1  2 
3  1  -1.01  1 
3  4-0.303135889  0  0 
subchain2  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4-0.303135889  0  0 
subchain3  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4-0.303135889  0  0 
subchain4  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4  -0.303135889  0  0 
subchain5  3  10 
3  1  -0.5  1  2 
3  1  -1.01  1 
3  4  -0.303135889  0  0 
branch6  2  1 
subchainl  310 
5  1  -0.5  1  2 
5  1-1.0  1  1 
5  4  -0.303135889  0  0 
subchain2  310 
5  1  -0.5  1  2 
5  1-1.0  1  1 
5  4  -0.303135889  0  0 
branch  7  21 
subchainl  3  10 
5  1  -0.5  1  2 
5  1-1.0  1  1 
5  4-0.303135889  0  0 
subchain2  3  10 
5  1  -0.5  1  2 
5  1-1.0  1  1 
5  4-0.303135889  0  0 
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Figure  3.  Bead-spring  dendrimer/comb  polymer.  Different  colors  correspond  to  the 
different  types  of  beads.  The  picture  corresponds  to  example  3.  Some  side 
chains  are  built  to  be  rigid. 


4.  Equilibration  of  the  Initial  Structures 


Short  simulation  runs  with  soft  potential  used  in  DPD  have  been  performed  for  an  ensemble  of 
polymer  chains  with  the  correct  end-to-end  distance.  DPD  potential  allows  chains  to  pass 
through  each  other  to  speed  up  the  polymer  dynamics.  Simulation  with  DPD  potential  deforms 
the  chains  on  the  intermediate  length  scales.  However,  this  deformation  is  fully  resolved  after 
500  tdpd-  After  this  step,  the  following  procedures  are  adopted  for  two  coarse-grained  models  1) 
Kremer-Grest  model  ( 4 )  and  2)  DPD  model  with  Segmental  Repulsive  Potential  (mSRP)  (5) 

Kremer-Grest  model.  Next  step  is  a  gradual  increase  of  the  strength  of  DPD  potential  in  order 
to  increase  distance  between  nonbonded  monomers  or  “push-off’.  Finally,  DPD  potential  is 
substituted  to  full  Lennard-Jones  potential.  This  step  does  not  significantly  alter  the  structural 
properties  of  the  melt.  Finally,  the  melt  is  simulated  with  Lennard-Jones  potentials  for  time  equal 
to  10000  tLj  to  remove  DPD  alterations  of  structural  properties  on  the  short  length  scale. 

DPD  with  mSRP  model.  mSRP  and  bending  potential  (see  section  2)  are  added  and  the  melt 
are  equilibrated  for  time  =  1000  topd/ihSrp  to  remove  DPD  alterations  of  structural  properties  on 
the  short  length  scale. 

To  equilibrate  the  initial  stricture  one  has  to  use  LAMMPS  software.  To  run  LAMMPS,  a  user 
has  to  copy  data,  chain  produce  by  gel  program  in  a  working  directory.  The  following  are  input 
scripts  for  LAMMPS  for  Kremer-Grest  and  DPD  with  mSRP  models,  respectively. 
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Kremer-Grest  model. 

units  Ij 
atomjstyle  bond 
special  bonds  Ij/coul  Oil 
read_data  data,  chain 
neighbor  0.4  bin 

neighjnodify  every  1  delay  1 
communicate  single  vel  yes 
bond_style  fene 
bond_coeff  1  30.0  1.5  1.0  1.0 

dump  mydump  all  custom  100000  chain _un  *.xyz  id  type  xu  yu  zu 
fix  1  all  nve 

timestep  0. 01 

thermo  1000 

pair_style  dpd  1.0  1.0  122347 

pair_coeff  *  *  25  4.5  1.0 

velocity  all  create  1.0  1  7786140 
run  50000 

write_restart  my_restart_dpd 
pair_coeff  *  *  50.0  4.5  1.0 
velocity  all  create  1.0  15086120 
run  50 

pair_coeff  *  *  100.0  4.5  1.0 
velocity  all  create  1.0  15786120 
run  50 

pair_coeff  *  *  150.0  4.5  1.0 
velocity  all  create  1.0  15486120 
run  50 

pair_coeff  *  *  200.0  4.5  1.0 
velocity  all  create  1.0  1  7986120 
run  100 

pair_coeff  *  *  250.0  4.5  1.0 
velocity  all  create  1.0  15006120 
run  100 

pair_coeff  *  *  500.0  4.5  1.0 
velocity  all  create  1.0  15087720 
run  100 

pair_coeff  *  *  1000.0  4.5  1.0 
velocity  all  create  1.0  15086189 
run  100 

write_restart  my_restart_dpdl 

pair_style  hybrid/overlay  Ij/cut  1.122462  dpd/tstat  1.0  1.0  1.122462  122347 

pairjnodify  shift  yes 

pair_coeff*  *  Ij/cut  1.0  1.0  1.122462 

pair_coeff  *  *  dpd/tstat  4.5  1.122462 

velocity  all  create  1.0  15086120 
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run 


50 

velocity  all  create  1.0  15086121 
run  50 

velocity  all  create  1.0  15086111 
run  50 

write_restart  my_restart _push 
velocity  all  create  1.0  15086125 
run  1000000 

write_restart  my_restart_lj 

DPD  with  mSRP  model. 

units  Ij 

atomjstyle  angle 
specialjbonds  Ij/coul  111 
newton  on  off 
read  data  data,  chain 
neighbor  2.5  bin 
neighjnodify  every  1  delay  1 
communicate  single  vel  yes 

#  syntax  for  bond-bond  pair_style 

#  use  hybrid/overlay  with  regular  dpd 

#  dpdbb  <cutoff>  <distance:  min=0  midpoint=2>  <search  for  top  volations:  no=0ye 
s=l>  <skin  distance> 

bond_style  harmonic 
bond_coeff  1225.0  0.85 
pair_style  dpd  1.0  1.0  3  73692 

pair_coeff  **  60.0  4.5  1.0 

velocity  all  create  1.0  1  7786140 

dump  mydump  all  custom  10000  chain _dpd*.xyz  id  type  xy  z 

thermo  1000 

fix  1  all  nve 

timestep  0. 01 

run  50000 

write_restart  my_restart_dpd 

velocity  all  create  1.0  40783012 

undump  mydump 

pair_style  hybrid/overlay  dpd  1.0  1.0  373692  dpdbb  0.8  2  0  2.5 

#  syntax  for  bond-bond  pair_coeff 

#  bondtype  bondtype  aO  cutoff  (optional) 

pair_coeJf  *  *  dpdbb  *  *  100.0  0.8 
pair_coeff  *  *  dpd  60.0  4.5  1.0 
angle_style  cosine 
angle _coeff  1  2.0 
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dump  mydumpl  all  custom  10000  chain*.xyz  id  type xy  z 
restart  10000  periodic _restart 
run  100000 

write_restart  my_restart_l 


5.  Results 


In  order  to  characterize  the  chain  conformations  in  the  melt  after  push-off,  the  authors  performed 
calculations  of  mean  square  internal  distances  of  chains  (R2{nfj/  n  is  averaged  over  all  possible 

combinations  of  segments  of  size  n  =  \i  -  j\  along  the  chains,  where  i  <  j  e  [1.JV]  are  monomer 

indices.  This  metric  is  an  excellent  indicator  of  chain  configuration  at  all  length  scales.  Code 
msid_l.c  is  used  to  calculate  this  function.  Arguments  are  N,  a,  Ntot,filel,file2,  where  A  is  a 
number  of  monomers  per  chain,  a  is  a  size  of  the  cubical  periodic  box,  Ntot  is  a  total  number  of 
particles,  filel  is  an  input  name  and  file2  is  an  output  name.  For  example,  ./msid  l  500  66.5 
250000  snap  lj.xyz  ljmsid.  Note:  input  files  should  be  in  VMD  xyz  format.  To  show  that  our 
preparation  technique  produced  a  sample  with  topology  equivalent  to  the  well-equilibrated 
structure,  we  applied  our  method  to  a  melt  of  M=  500  chains  of  N=  500  beads.  Simulation  box 
was  66.5  x66.5><66.5  a  .  We  compared  a  topology  of  1)  the  melt  produced  by  our  method  or  2) 
the  melt  equilibrated  by  running  MD  simulation  for  time  long  enough  that  the  chains  moves  3 
times  their  own  size  (the  run  time  was  t  =  6.5  x  106  zu  ).  Figure  4  shows  the  results 
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Figure  4.  Evolution  of  mean  square  internal  distances  for  a 

melt  of  chain  length  N  =  500  throughout  our  preparation 
procedure.  Blue,  red  and  green  symbols  correspond  to  the 
chains  after  initial  1 000  xDPD  simulation  with  the  soft  (DPD) 
potential,  after  fast  “push-off’  and  after  1 0000  xLJ  simulation 
with  Lennard-Jones  potential,  respectively.  The  black 
symbols  denote  our  target  function  obtained  for  the  equilibrated 
melt,  where  simulation  time  is  t  =  6x106tlj  . 

Using  geometrical  analysis,  the  authors  estimated  the  entanglement  length,  Ne  by  running 
Zlcode  (9).  The  authors  found  that  Ne  for  an  equilibrated  and  prepared  system  is  83  and  84, 
respectively. 

The  results  demonstrate  the  excellent  agreement  between  structural  properties  of  the  melts 
generated  via  our  method  and  brute-force  equilibration. 
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